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DETERMINATION OF THE MASS LOSS RATE AND THE TERMINAL 
VELOCITY OF STELLAR WINDS. I GENETIC ALGORITHM FOR 
AUTOMATIC LINE PROFILE FITTING 

Leonid Gcorgicv, 1 Xavicr Hernandez, 1 
RESUMEN 

Se presenta un nuevo metodo de ajuste automatico de perfiles P Cyg observados en espectros UV de vientos 
estelares. La funcion fuente de las hneas se calcula usando la aproximacion de Sobolev y el flujo emitido se 
obtiene atraves de la solution exacta de la equation de transporte (similar al metodo SEI descrito por Lamer 
et al. 1987). La calidad del ajuste se evalua usando el cstimador Likclyhood. La maximisation del Likelyhood 
se hace atraves de un algoritmo genctico. Las ventajas de nuestro metodo comparado con otros metodos 
similares son su cstabilidad y su poca sencibilidad de las condiciones iniciales. Ademas, el metodo garantiza 
la localization del maximo global de la superficie del Likelyhood, cual no se puede hacer con los algoritmos 
clasicos. Presentamos la implementation del metodo, las pruebas de su funcionalidad usando datos sinteticos 
y reales y la estimation de los rangos de confianza de los resultados. 

ABSTRACT 

A new method for automatic fitting of P Cyg line profiles in UV spectra of stellar winds is presented. The 
line source function is calculated using Sobolev approximation and the emergent flux is obtained by exact 
integration of the equation of the radiation transport (similar to the SEI method described by Lamers et al. 
(1987)). The quality of the fit is evaluated using the Likclyhood estimator. The maximization of the Likelyhood 
is done by a genetic algorithm. The advantages of our method with respect to other similar approaches are its 
robustness and its insensibility to the initial guess. In addition, the algorithm guarantees the localization of the 
global maximum of the Likclyhood hypersurface, which is not the case for classical minimization algorithms. 
Here we present an implementation of the genetic algorithm for line profile fitting, its tests on both synthetic 
and real data and and estimation of the confidence limits of the results. 
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1. INTRODUCTION 

Many hot stars lose mass through a supersonic 
wind. The terminal velocity of these winds, and the 
details of the mass loss rate play an important role 
in stellar evolution. Stellar winds crucially determine 
the interaction of the star with its surrounding inter- 
stellar medium, this being one of the central feedback 
mechanisms, which in turn modulate star formation 
at galactic levels. Estimating stellar mass loss rates 
is therefore a highly relevant problem. 

The best way of determining the mass loss rate 
of a star is by studying the radio emission from 
the wind itself (Lamers et al. 1999) however, with 
current observational techniques this is limited to 
nearby stars with high mass loss rates. The other re- 
liable method is based on a detailed modeling of the 
stellar spectrum. The recent development of realis- 
tic codes which include many atomic lines makes the 
latter approach the preferred one. Unfortunately, 
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the fitting of the spectrum requires high quality 
data which covers a wide range of wavelengths and 
with high dispersion, this again limits the method 
to bright or nearby stars. An alternative has been 
proposed by Lamers et al. (1987) and Grocnewc- 
gen & Lamers (1999). This method is based on a 
code which approximates the opacity of the wind 
as a simple function of the distance from the stel- 
lar nucleus. The code calculates the source func- 
tion in a given line within the Sobolev approxima- 
tion (Sobolev, 1957) and then calculates accurately 
the emitted profile. The fit of the observed profile 
gives the terminal velocity of the wind together with 
the total optical depth along the line of sight. The 
mass loss rate is estimated from the later (see next 
section). Even though the obtained value is sensitive 
to certain assumptions regarding the ionization frac- 
tion of the element whose line is fitted, the method 
has the advantage of being applicable to faint stars 
for which only low resolution and low signal to noise 
spectra are available. 
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The main disadvantage of this approach is the 
large number of free parameters which have to be 
adjusted. Groenewegen & Lamers (1989) proposed 
a method for determination of the parameters by a 
"classical" minimization of the \ 2 = (1 — F Q /F C ) 
where F a is the observed flux and F c is the calcu- 
lated one. But the multi-dimensional fitting con- 
verges to the correct solution only if the initial guess 
is close to the \ 2 minimum. The complex topology 
associated with the high dimensional hypersurfacc 
being treated leads to the appearance of numerous 
local minima. The "classical" minimization meth- 
ods usually get trapped in the local minimum close 
to the initial guess missing the global one. The some- 
what lengthy numerical procedure involved, again 
together with the large number of dimension of the 
problem, makes a dense sampling of the parame- 
ter space completely impractical. Recent years have 
seen the development of non-standard maximization 
techniques such as simulated annealing and genetic 
algorithms, particularly in connection with astro- 
physical applications (e.g. Teriaca et. al 1999 in 
fitting radiative transfer models to solar atmosphere 
data; Sevenster et al. 1999 in fitting parameters of 
Galactic structure models to observations or Met- 
calfe 2003; for the fitting of White dwarf astroseis- 
mology parmeters) . These methods provide a frame- 
work which in a reliable manner identifies the global 
maxima (or minima, as the case might be) of com- 
plex multi-dimensional surfaces, with a close to op- 
timal use of computing resources. 

In this paper we present an application of a ge- 
netic algorithm used to obtain the optimal param- 
eters of a stellar wind, within the Sobolev approxi- 
mation, for the fitting of P Cyg line profiles of UV 
resonance lines. 

Section (2) describes the radiative transfer 
model. In section (3) a description of the genetic 
algorithm is given, with use of it in a number of syn- 
thetic test cases and real data examples appearing 
in section (4). Section (5) presents our conclusions. 

2. THE ALGORITHM 
2.1. Radiative transfer code 
The core of the method is a code which solves the 
radiative transport problem under the assumption of 
a spherically symmetric stellar wind. Following the 
SEI code (Lamers et al. 1987) we approximate the 
radial Sobolev optical depth t(v) as: 
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¥Q is the velocity of the wind, nor- 



opacity of the line, vq is the laboratory frequency of 
the transition and c is the speed of light. T, ct\ and 
cti are free parameters of the model, which describe 
the radial dependence of the optical depth. The ve- 
locity field of the wind is adopted to be a standard 
/3-law 
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where is the terminal velocity. The winds 
of the hot stars are not smooth but have random 
motions. The shape of the blue wing of the P Cyg 
absorption component suggest that the intrinsic line 
profile has a width on the order of hundreds of kilo- 
meters per second. The physical nature of these 
random motions is not clear but they act as an ad- 
ditional line broadening mechanism similar to tur- 
bulence and are frequently called "turbulence" al- 
though they may have little to do with it. We model 
this chaotic motions assuming the intrinsic line pro- 
file to be Gaussian defined as 



(v) = 



1 



-exp 



(3) 



rplVturb 

The turbulent velocity, Vturb, can be variable 
throughout the wind, but because its physical na- 
ture is not clear there is no model of its possible 
variability. To keep the model as simple as possible, 
we assume Vturb to be constant and a parameter of 
the model. 

Finally, we set the innermost point of the grid, 
Rmin to the location where V(R m i n ) is equal to the 
sound speed V soun d- This is done by setting the pa- 
rameter Rq in J2J) as 



Ro — Rn 
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v sound * 
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Equations Q, J5J and Q determine the distri- 
bution of the central opacity of the line Xo as a func- 
tion of the radius r. Within the Sobolev approxima- 
tion, one can now calculate the line source function 
S and then can apply a formal solution of the equa- 
tion of radiative transport throughout the wind and 
calculate the emergent flux (see Georgiev & Kocnigs- 
berger, 2004 for more details of the code). The code 
is designed to work in 3D geometry, but for the cur- 
rent test purposes, the solution is restricted to the 
case of spherical symmetry. Finally, we calculate the 
mass loss rate as follows (Groenewegen & Lamers, 
1989). The line opacity xo is 



malized to the terminal velocity , xo is the central 



ire- 

Xo = Jiv.ni 



t(v) 



vq dV(r) 
c dr 
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where n; is the population of the lower level, fi u 
is the oscillator strength and we do not take into 
account the correction for stimulated emission. The 
level population m can be written as 



Ni on q ex = N atom q ion q ex — NnA E q. 



tonnex , 



(0) 

where qi on and q ex are the ionization and excitation 
fractions and Ae is the chemical composition of the 
element relative to Hydrogen by number. Then the 
level population ri\ is related to the mass loss rate 
by the equation of continuity 



ni 



A E Mq ion q ex 
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where /i is the average molecular weight. Substitut- 
ing Q into JSJl, one can obtain the product Mq ion q ex 
at each point r as 



Mq t , 



dV(r) a 
— ; — r V(r) 
dr 
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Then, if the parameters I4o, Vt ur b, P, T, ct\ and 
«2 are given, equation and © together with the 
atomic data and the chemical composition allow the 
estimation of the product Mqi on q ex . The method 
described in this paper suggests a better way to ob- 
tain the six parameters that describe the velocity law 
and the optical depth. It may not improve the preci- 
sion of the calculated mass loss rate if the dominant 
uncertainty is the ionization fraction. 

Applying the above procedure, we are able to 
calculate any line profile, once the six parameters 
of the method V^,, Vturb, P, T, a\ and cti have been 
specified. The next step is the selection of a statisti- 
cal estimator of the goodness of fit to be associated 
to a particular set of parameters and a given ob- 
served profile. From a Bayesian perspective, the op- 
timal such indicator is the Likelyhood, calculated as 
the probability that a certain data set be observed, 
given a model. We hence compare a modeled pro- 
file against a given observation through a Likelihood 
estimator defined as: 



logL(V 00 ,V t urb,l3,T,a 1 ,a2)= (9) 
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Where F i o and F i c are the observed and calcu- 
lated flux respectively at a certain wavelength, \ 
and the factors Ui are the observational errors at the 
same A;. The weight a assigned to each observed 



point is taken from the signal-to-noise ratio of the 
observations. Setting a large Oi to certain points 
forces them to participate less in the final fit, which 
allows us to exclude effectively spectral regions where 
overlap with other lines occur. Hence, the final re- 
sult depends on the assigned er^. In practice, the 
bands of the observed spectrum affected by spurious 
details are easily detectable and the assignment of 
the weights is straightforward. 

The synthetic profile is calculated over a grid ad- 
justed for better sampling the velocity field in the 
wind. The calculated profile is then interpolated on 
the wavelengths of the observed profile using mono- 
tonic cubic polynomials (Steffen, 1990). Finally, 
equation © is calculated and used as a measure of 
the quality of the fit. 

2.2. The genetic algorithm 

Once we have a method for calculating the line 
profile which corresponds to a particular set of six 
parameters, and an estimator of the goodness of the 
fit, we must now determine what the optimal param- 
eter vector is; i.e. we must find the combination of 
six parameters which maximizes our Likelihood func- 
tion. This will yield the physical parameters of the 
stellar wind associated to the observed line profile, 
from which the mass loss rates and wind terminal 
velocities can be estimated. 

A dense exploration of our six-dimensional pa- 
rameter space is unfeasible. Classical maximization 
techniques are also unreliable in connection to prob- 
lems of such high dimensionality, and liable to get 
trapped by local maxima (Press et al . 1989). In 
fact we have encountered numerous test cases where 
local maxima exist in the Likelyhood surface of our 
problem. We bypass these obstacles using genetic al- 
gorithms in searching for the maximum of the Like- 
lyhood surface. 

The idea is to simulate a population of organisms 
(we call them " bugs" ) which breed and evolve follow- 
ing prescriptions based on that biological systems are 
thought to follow, the result being a progressive in- 
crease in the fitness of the population, with the fittest 
individual in each generation eventually reaching the 
absolute maximum of the fitness surface. For the 
case at hand, a "bug" is a set of wind parameters 
described above. 

Once an observed line profile has been picked, 
the first step is to select TV random points in our pa- 
rameter space, each a parameter vector for the stel- 
lar wind, (Voo, Vturb, P, T,a\, a2). The ranges over 
which these parameters are to be chosen have to be 
determined by inspection of the observed line profile. 
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Only ample margins enclosing the expected value are 
needed. We turn to this issue in the examples sec- 
tion. 

The six coordinates of each point are then turned 
into binary numbers, with a pre-determined resolu- 
tion and then combined in a string. We have used 8 
bits per parameter for our implementation, resulting 
in a string of Nq = 6 parameters x8 bits = 48 bits. 
This selection gives to our search algorithm a res- 
olution of Rangei/2 8 in each of the six dimensions 
where Rangeis are the ranges for each parameter 
chosen above. 

Each of the TV randomly selected points becomes 
a "bug" of the first generation. The string of TVq I s 
and s which corresponds to its 6 coordinates in pa- 
rameter space becomes its "genotype" , and the Like- 
lyhood associated to that point becomes its "pheno- 
type" . The TV members of the first generation are 
then ranked according to the value of the Likelihood 
associated to each, with those ranking above N/Sp 
being deemed "superior" , and those ranking below 
this threshold "inferior". 5p is a parameter of the 
simulation which has to be a number greater than 
1.0, typically around 1.5. 

This completes the first generation, which arc 
now "mated" to produce the second one. This pro- 
ceeds by the random selection of pairs of individuals 
from the pool of TV members of the first generation. 
Once a pair is selected, "offspring" are calculated. 
The first Cp genes of one of the parents are taken 
for first Cp genes of the offspring, and the remain- 
der are taken from the other parent. Cp is a "cross- 
ing point", chosen at random in the interval (0-TVg), 
each time an offspring is constructed. 

In this way, the offspring will have "genes" in 
some aspects resembling those of one parent, and in 
some cases the other, the corresponding position in 
parameter space will hence share some coordinates 
with one parent, and some with the other, with one 
of them being a mixture. Also a fraction Mp of 
randomly selected genes in the new generation are 
"mutated" i.e. they are flipped from 1 to 0, or to 
1, as the case might be. Mp is a second and last 
parameter of the method, which together with Sp 
must be chosen and tuned through some testing of 
the algorithm in the context of the particular prob- 
lem being treated. We have used Mp close to 0.01 
which means 1% of the genes are mutated 

The number of offspring to be produced by each 
pair depends on the fitness of the parents, if both 
are "superior" individuals, 3 offspring are calculated. 
The pairing of a "superior" and an "inferior" individ- 
ual yields only 2 offspring, and when two unfortunate 



"inferiors" mate, only one child ensues. The process 
is repeated until N offspring have been obtained, the 
new generation has been constructed, and it com- 
pletely replaces the former. 

The parameter Sp hence represents a selection 
pressure. If chosen close to 1, most individuals will 
be deemed "superior" , and hence the second gener- 
ation will be distributed in parameter space much 
like the first one was. However, as Sp takes larger 
values, the second generation will be made up of the 
sons of (mostly) the "superior" individuals of the 
former, yielding a gradual climb up the Likelihood 
surface of the problem. Next the genotypes of the 
new generation are turned back into coordinates in 
our parameter space, and then into phenotypes -the 
associated Likelihood. The ranking is repeated, and 
the cycle iterated. 

In this way, we have introduced mating, selec- 
tion, and chance mutations, the three main ingredi- 
ents of biological evolution, albeit in a highly sim- 
plified manner. Selection forces subsequent gener- 
ations up the Likelihood surface, with the random 
mutations ensuring that a fraction of the individuals 
are constantly sampling new regions, which in turn 
guarantees the absolute maximum will eventually be 
found. 

In many classical maximization algorithms one 
is called to evaluate the gradient of the surface at 
any given point. When the surface is not an analyti- 
cal one, but the result of a lengthy numerical proce- 
dure, as in our case, the gradients could quite easily 
be dominated by roundoff errors. The method de- 
scribed calls for no gradients of the surface, only its 
value at the points being tested. With the resolution 
described here, a dense sampling of our Likelihood 
surface would have required (2 8 ) 6 = 2.8 x 10 eval- 
uations of our wind model. We have used N=100 
members per generation, and both in synthetic ex- 
amples where the answer was known in advance, and 
in tests with real data, have found convergence in 
200 generations, i.e. only 2 x 10 4 costly line profile 
calculations. 

The method described here differs only slightly 
from the general purpose genetic maximization al- 
gorithm PIKAIA, described in Charbonncau (1995), 
and having been used successfully in a number of as- 
trophysical applications to date (e.g. Teriaca et. al 
1999, Metcalfe 2003). The differences between the 
two are limited to the way "selection" is treated, the 
implementation described here giving results more 
suitable for the particular problem we are treating. 
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Fig. 2. The ellipses show 1 — a confidence intervals on the 
recovered parameters (circles) for the fit to the synthetic 
profile with S/N=100. The input values are shown by 
the solid dots. 

3. TEST CASES 

The first test of the algorithm is made on syn- 
thetic data. We calculate the line profile of the 
C IV 1548/1550A doublet with a predefined velocity 
law and wind density distribution and then fit this 
profile using the genetic algorithm. Fig. 2] shows the 
original profile and the resulting fit for three differ- 
ent signal-to-noise ratios. The parameters used to 
generate the profile and the result of the fitting are 
summarized in Tab.^ 

The genetic algorithms do not include an easy 
way of estimating the errors of the results. The 
method only serves to locate the global maximum in 
an efficient manner. We estimate the confidence in- 
tervals of the recovered parameters by Monte-Carlo 
simulations of synthetic line profiles. In this way, all 
uncertainties inherent to the method and the mod- 
elled data are taken into account. We select three 
signal-to-noise ratios: 100, 50 and 10. For each of 
them we generate 20 different line profiles, using the 
same synthetic data but adding different and inde- 
pendent noise. The obtained line profiles were fitted 
by the genetic algorithm using different first genera- 
tions of "bugs" . Figs . 12131 and EH show the confidence 
intervals obtained through this method for the recov- 
ered parameters for the simulated profiles of Fig. ^ 

Once the Monte-Carlo simulations are all done, 
we find the means for all the recovered parameters, 
and of the recovered mass loss rates. Next, a full co- 



Fig. 3. The ellipses show 1 — a confidence intervals on the 
recovered parameters (circles) for the fit to the synthetic 
profile with S /N=50. The input values are shown by the 
solid dusts. 
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Fig. 4. The ellipses show 1 — a confidence intervals on the 
recovered parameters (circles) for the fit to the synthetic 
profile with S/N=10. The input values are shown by the 
solid dusts. 

variance analysis is performed, yielding the param- 
eter's standard deviations, which together with the 
relevant covariances can be used to construct statis- 
tically meaningful confidence ellipses. Fig. [5] shows 
the 1 — a ellipses for the case of S/N =100, with the 
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Fig. 1. Test of the fitting algorithm. A synthetic profile of C IV 1548/1550A (squares) is shown together with the best 
fit (solid line) for three adopted S/N ratios. 



solid dots indicating the input values, and the circles 
showing the ones reached by the fitting algorithm. 
The mass loss rate is recovered almost exactly, with 
a 1 — a interval of 30%. A negative correlation is 
seen between the inferred M and the recovered ter- 
minal wind velocity, while mass loss rate and (3 show 
a weak positive correlation. Figs. Eland 01 are equiva- 
lent to Fig. |3 but for the cases with S/N values of 50 
and 10, respectively. The same qualitative features 
can be seen, with the ellipses getting only slightly 
larger, to reach 40% for the S/N =10 case. It is re- 
assuring that errors on the inferred mass loss rate are 
not highly sensitive to the value of the S/N of the 
simulated profile, making the method usefull even in 
cases where observations are not of the best quality. 
Also, we note that no systematics appear on the re- 
covered M, we have constructed unbiased estimator 
of the important physical parameter we set out to in- 
fer. We have to stress that the error on the mass loss 
rate are the formal errors of the fit. The errors re- 
lated to the determination of the ionization fraction 
will be treated in the next paper of the series. 

As can be seen from Tab. ^ m all cases the al- 
gorithm finds the correct maximum of the likelihood 
estimator, the mass loss rate is recovered success- 
fully, always well within 1 — a of the input value, 
even at low S/N values. As an additional check we 
tried to fit the line profiles of real stars. Fig.|Slshows 
the profile of C IV 1548/1550A doublet observed in 
IUE spectrum of £ Pup (HD668110) and Fig.dshows 




-4000 -2000 2000 4000 

Velocity [km/a] 



Fig. 5. Fit of the C IV 1548/1550A line in ( Pup 
(=HD66811). The thin line is the observed spectrum. 
The dotted line is the fit 

the profile of Si IV 1398/1402A doublet in the IUE 
spectrum of HD 30614. The comparison between 
the results of our fit and the parameters obtained 
by Groenewegen & Lamers (1989) (Tab. GJl shows a 
good agreement. Our solution shows a slower veloc- 
ity law (large (3) but our final fit has higher Likely- 
hood than the fit obtained with the published pa- 
rametes. For both stars we obtain a larger and 
lower Vturb- As shown in Figs . [HI the two parameters 
are correlated and one could expect this behavior. 
But in both cases the automatic fit gives values very 
close to the more human controlled solution. 
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TABLE 1 



PARAMETERS OF THE SYNTHETIC PROFILE AND THE RESULTS OF THE FITS. 



Parameter 


Exact value 


S/N=100 


S/N=50 


S/N=10 


V x (km s" 1 ) 


2000 


1997 ± 19 


2007 ± 23 


2014 ± 30 


Vt U rb (km s _1 ) 


100 


114 ± 12 


107 ± 13m 


101 ± 9 


(3 


1.0 


0.99 ± 0.06 


0.99 ± 0.05 


0.87 ± 0.08 


T 


2.0 


2.01 ± 0.06 


1.96 ± 0.08 


2.015 ± 0.11 


log Mq ion q ex (M G /yr) 


-7.0 


-7.00 ± 0.13 


-7.00 ± 0.14 


-7.02 ± 0.15 



TABLE 2 

PARAMETERS OF THE WIND OF HD30814 AND HD66811 DETERMINED IN THIS WORK 
TOGETHER WITH THE SAME DATA OBTAINED BY GROENEWEGEN & LAMERS (1989). 



Parameter 


HD66811 C IV 1548/1550A 


HD30814 Si IV 1398/1402A 




This work 


Published 


This work 


Published 


Voo (km s- 1 ) 


1550 ± 16 


1550 ± 50 


2137 ± 21 


2200 ± 60 


Vturb (km s" 1 ) 


240 ± 24 


190 ± 70 


287 ± 30 


290 ± 70 




1.13 ± 0.1 


0.7 ± 0.1 


1.18 ± 0.1 


0.7 ± 0.1 


log Mq lon q e x (M Q /yr) 


-6.99 ± 0.15 


> -7.14 


-7.27 ± 0.15 


> -8.4 
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Fig. 6. Correlation between V x 
synthetic profile with S/N=50. 



and Vturb calculated for 



4. CONCLUSIONS AND DISCUSSION 

The determination of the basic stellar parame- 
ters like wind velocity and mass loss rate is a very 
time consuming process. In this paper we present a 
method which allows an objective determination of 
these parameters based on a completely automatic 
procedure. The test cases shown above support the 
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Fig. 7. Fit of the Si IV 1392/1402A line in HD 30614. 
The dotted line is the fit 



robustness of the method even in case of low S/N 
data. The method is especially useful for objects 
for which only low resolution and low S/N data are 
available. Their spectra show only few spectral fea- 
tures and the application of full NLTE codes is not 
practical. Our method was successfully applied by 
Arrieta & Stangellini (2004) for central stars of LMC 
planetary nebulae. Also, the use of Monte-Carlo 
simulations on the recovered parameters allows for 
the secure and objective determination of statisti- 
cally meaningfull confidence intervals around recov- 
ered parameters, indeed, the full covariance matrix 
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Fig. 8. Fit of the P V 1117/1129 line in the LMC star 
HD5980. The thin line is the observed spectrum. The 
dashed line is the fit. Small weight were applied to the 
interstellar absorption line, which contaminate the P Cyg 
profile. 

is available. But. we have to stress, that the parame- 
ter obtained by this method is not the mass loss rate 
M itself, but rather the product Mq ion q ex . There 
is no easy way to separate these quantities on the 
base of low dispersion and low S/N data on which 
only few lines are measurable. Our method helps to 
obtain a more reliable estimates of the product in an 
easier way, but the final decision on how to disen- 
tangle M and the ionization and excitation fractions 
depends on the studied object and on the problem 
which needs to be solved. 

An important advantage of our method, com- 
pared to other fitting algorithms is a minimal effect 
of the initial guess on the final solution. This is be- 
cause the method does not rely on a series expansion. 
The genetic algorithm requires only a crude margins 
on the parameters to find the maximum of the likc- 
lyhood surface. In case the maximum is not within 
the proposed ranges, the population migrates to the 
parameter's limit which is closer to the maximum. 
A subsequent run with new ranges that increase this 
limit usually leads to a successfully solution. The 
price to pay is only an increased number of pro- 
file evaluations. But the fast computers make this 
price affordable. The calculation of one generation 
of "bugs" takes less than 2 minutes on a Pentium IV 
class computer. Usually fifty to a hundred genera- 
tions are needed to fit one line profile. Each "bug" in 
a generation is independent of the rest of the popu- 
lation, so a simple parallelization increases the speed 
of the calculation almost proportional to the number 
of available processors. 

The method has however few disadvantages. 
First of all, the choice of the weights Ui in equa- 
tion l^J is not automatic. The observed profiles fre- 



quently have overlapped absorption lines of other el- 
ements, interstellar lines or absorption components 
of the same star and the same element but which are 
not formed in the wind. The weights <ii should be 
chosen lower in the affected bands so the Likelyhood 
is calculated using only the pixels which belongs to 
the line. FigEJshow the fit of the P V 1117/1128A 
line in the star HD5980 as an example of the per- 
formance of the code on a highly contaminated line. 
The FUSE spectra are characterized by a huge num- 
ber of interstellar absoption line. In the above case, 
a small weight (large errors) were assigned to the 
bands with high contamination. As seen from FigEl 
the code finds a satisfactory fit which reproduce the 
P Cyg line profile, avoiding the contamination. 

Another problem is that the parameters can be 
obtained only from a grid. As a result, the code never 
converges to the exact solution but to the node of the 
parameter's grid which is closest, but not exactly at, 
the maximum of the likelyhood surface. That might 
cause the systematic error seen in Fig. 0] One pos- 
sible solution to this problem is a second run of the 
code with smaller ranges of the parameters around 
the solution found. Another approach is to com- 
bine a genetic algorithm with a classical minimiza- 
tion method. The genetic algorithm guarantees, that 
the solution is near the global maximum of the like- 
lyhood surface and then the classical minimization 
algorithm starts from that point and finds the exact 
solution. Until now we have explored only the first 
option. The combination of the algorithms and their 
application to a wide range of objects, together with 
the problem of separating the mass loss rate from the 
ionization fraction will be subject of the next paper 
of this series. 

We thank Gloria Koenigsberger for many use- 
full discutions and for critical reading of the 
manuscript. This work was supported by CONA- 
CyT grants 40864 and 42809 and UNAM/DGAPA 
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